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When a Bose-Einstein condensed cloud of atoms is given some angular momentum, it forms 
vortices arranged in structures with a discrete rotational symmetry. For these vortex states, the 
Hilbert space of the exact solution separates into a “primary” space related to the mean-field Gross- 
Pitaevskii solution and a “complementary” space including the corrections beyond mean-field. Con¬ 
sidering a weakly-interacting Bose-Einstein condensate of harmonically-trapped atoms, we demon¬ 
strate how this separation can be used to close the conceptual gap between exact solutions for 
systems with only a few atoms and the thermodynamic limit for which the mean-field is the correct 
leading-order approximation. Although we illustrate this approach for the case of weak interactions, 
it is expected to be more generally valid. 

PACS numbers: 03.75.Lm, 05.30.Jp, 67.85.Hj, 67.85.Jk 


I. INTRODUCTION 

Cold atomic quantum gases are typically dilute with 
an average interatomic distance much larger than the 
scattering length for atom-atom elastic collisions. This 
justifies the use of the mean-field approximation, which 
assumes a simple product form for the many-body wave- 
function in the case of bosonic atoms. The complicated 
many-body problem is then reduced to one of a single 
variable, and the effect of interactions is described by a 
nonlinear term. This procedure was developed by Gross 
and Pitaevskii some decades ago [l|, 0 . Since the experi¬ 
mental realization of Bose-Einstein condensation (BEG) 
in trapped atomic gases, this approach has been used 
with remarkable success, see e.g., Refs. [^-Q. 

One of the many fascinating effects associated with 
the superfluid properties of these gases is the formation 
of vortices in response to rotation. When the ratio of an¬ 
gular momentum to particle number increases, the num¬ 
ber of vortices in the cloud grows, and they group in 
structures with discrete rotational symmetries (as illus¬ 
trated by the mean-held densities in Fig. [T] below). Such 
vortex states have been observed in a number of experi¬ 
ments, see for example. Refs. The literature on 

this topic is extensive, as summarized by the reviews 
(isl - fl^ . For a dilute and harmonically-trapped Bose- 
Einstein condensate of atoms, the rotational properties 
have been thoroughly analyzed both within the Gross- 
Pitaevskii approximation as in Refs. [ira , and bey ond 
the mean-held approximation as in Refs. 1^ |2(R^ . 

Going beyond the mean-held approximation one often 
applies the so-called conhguration-interaction (Cl) for¬ 
malism. In this numerical approach, one typically uses 
the Fock states constructed from a given set of single¬ 
particle states as a basis for the expansion of the exact 
many-body wavefunction. Other approaches are often 
variational, such as quantum Monte-Carlo or density 
functional techniques for correlated Bose gases . The 
so-called coupled-cluster approach, originally formulated 


for nuclei [13,13 often applied to atomic and molec¬ 
ular systems of fermions [13, S] has also been adapted 
to bosonic systems [5l| and is based on a series expan¬ 
sion of excitation operators acting on the corresponding 
mean-held ground state conhguration. 

An important advantage of the full Cl approach is that 
apart from an almost always inevitable truncation of the 
Hilbert space, no further assumptions are made regard¬ 
ing the functional form of the many-body wavefunction. 
The method fully accounts for the correlations between 
the particles and accurately describes the low-lying ex¬ 
citations. However, the dimension of the Hamiltonian 
matrix grows very rapidly with the number of particles. 
Thus, little is currently known about the intermediate 
regime between small and large systems, in which exact 
diagonalization becomes prohibitively difficult but the 
mean-held approach still suffers from signihcant hnite- 
size corrections. 

Here, we wish to shed new light on this problem, sug¬ 
gesting a procedure that offers direct insight into the 
question of how a hnite-size system of bosonic particles 
approaches the thermodynamic limit in which the Gross- 
Pitaevskii approach is known to be exact [H, 

With increasing particle number, we hnd a power-law 
convergence of the exact ground state into the mean- 
held Gross-Pitaevskii solution. Our study thus provides 
a clear and general strategy for this problem and offers 
strong arguments for both its validity and practicality. 
Further, these arguments are not limited by either the 
particle number or by the strength of the interaction. 


II. ROTATIONAL PROPERTIES WITHIN THE 
GROSS-PITAEVSKII APPROXIMATION 

We begin by briehy reviewing some of the well-known 
results regarding the formation of vortices in a weakly- 
interacting dilute Bose-Einstein condensate. The crite¬ 
rion adopted here for the formation of vortices is their 



2 


O ^ ^ 

FIG. 1: Iso-surfaces of the density distribution of the bosonic 
cloud within the mean-field approximation (taken at about 
one third of the maximum value) for different ratios of angular 
momentum to particle number, L/N = 1,1.8, and 2.3 from 
left to right (first discussed by Butts and Rokshar 0). The 
vortices appear as holes, and the phase of the order parameter 
(in colorscale from red to blue) jumps by 2n when encircling 
each vortex. 


energetic stability, namely the minimization of the en¬ 
ergy either for a fixed value of the angular momentum or 
for a fixed value of the rotational frequency of the trap. 
We consider a harmonic trapping potential that is very 
tight along the axis of rotation, here chosen to be the 
z axis, with oscillator frequencies oj = ujx = i^y ^ oJz. 
The two-body potential is assumed to be the usual con¬ 
tact interaction. For sufficiently weak interactions, fiWz 
is much larger than the interaction energy, and the atoms 
occupy only the ground state of the harmonic potential 
in the z-direction. Thus, the problem becomes effectively 
two-dimensional. The interatomic potential has the form 
F(ri —Tj) = itQ (5(ri —Tj), with uq = [/q J\(j){z)\‘^ dz. Here 
Uq = ^-nh^a/M is the matrix element for elastic colli¬ 
sions, a is the corresponding scattering length, M is the 
mass of the particles, and (j){z) is the ground state of the 
oscillator potential along the z axis. The Hamiltonian 
that we consider is thus 


H 


^ -t 2 \-7 2 T\ ,f ^ 

E Ti V,- J\d 2/2 2\ ^0 r/ \ 

1=1 i^j=l 


( 1 ) 


Within the mean-field Gross-Pitaevskii approxima¬ 
tion the many-body wavefunction is assumed to have 
a product form, while the corresponding order param¬ 
eter y) can be expanded in the eigenstates of the 

non-interacting problem. Making the assumption that 
the interaction energy is much smaller than hui, it is 
sufficient to consider only the single-particle states (j)o,m 
of the lowest Landau level with zero radial nodes and 
angular momentum mh > 0. (As we will explain in 
Sec. IV, the assumption of weak interactions is not es¬ 
sential, and the approach presented below is expected to 
remain valid for stronger interactions). The order pa- 
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C3 

C4 C6 

1.0 

0 

1.0000 

0 

0 

0 0 

1.8 

0.4803 

0 

0.7992 

0 

-0.3611 -0.0127 

2.3 

0.5312 

0 

0 

0.8179 

0 -0.2210 


TABLE I: Coefficients Cm of the expansion of the Gross- 
Pitaevskii order parameter in the lowest Landau level for 
0 < m < 6 (see text). 


rameter 'I'mf is thus expanded in the basis of the states 
4 > 0 ,m, ^'mf = where the amplitudes Cm 

are variational parameters. 

Minimizing the energy functional subject to the con¬ 
straint of fixed L/N (where Lh is the total angular mo¬ 
mentum, and N the number of particles in the trap), 
one finds that as L/N increases, there is a sequence of 
phase transitions associated with the formation of one 
or more vortices in the gas (see the work by Butts and 
Rokshar [2^ and, for example. Refs. [2ll - l^ lSj| ). This is 
a direct consequence of the fact that for a hxed value of 
L/N only certain single-particle states are occupied by 
a macroscopic number of atoms of 0{N). For example, 
for L/N = 1 the mean-field approximation yields a so¬ 
lution in which the only occupied single-particle state is 
the one with m = 1 and there is a single vortex at the 
trap center (see Refs. [^.[^. [261 - 1^ ). A similar behavior 
is found at higher values of L/N [^, [2^. For exam¬ 
ple, for L/N = 1.8, only the single-particle states with 
TO = 0, 2,4,... are occupied and the order parameter has 
two-fold symmetry. Correspondingly, for L/N = 2.3, the 
mean-held state consists of single-particle orbitals with 
only TO = 0, 3, 6 ,... and it has three-fold symmetry, 

The actual values of the variational parameters Cm 
which are derived within this method are given in the Ta¬ 
ble below for L/N = 1.0,1.8 and 2.3 (with 0 < to < 6). 
In Fig.[T] we show the density isosurfaces for these states 
(having one-fold, two-fold, and three-fold symmetry), 
where the color scale (from red to blue) indicates that 
the phase of the order parameter changes by 27r when 
encircling the vortex singularity . 


III. EXACT SOLUTIONS COMPARED TO 
MEAN-FIELD 


In order to obtain the exact solution of the prob¬ 
lem m, in, m ) one diagonalizes the many-body Hamil¬ 
tonian H subject to the constraints of fixed particle num- 
ber, fixed total angular momentum, 

Sm With increasing values of N and L, the 

exponentially growing computational complexity severely 
restricts the size of numerically tractable systems to a 
few dozen atoms at most. Yet, studying the detailed 
structure of the exact wavefunction suggests a substantial 
simplification: In analogy with the mean-field approach 
discussed above, when an yrast state (i.e., the state with 
lowest energy at fixed angular momentum L) has a given 
discrete rotational symmetry, only certain single-particle 
states are occupied by a macroscopic number of atoms 
of 0{N). This permits a separation of the total Hilbert 
space into a “primary” subspace that includes only Fock 
states constructed exclusively from single-particle states 
which are macroscopically occupied within the mean-field 
approximation and a far larger “complementary” space 
that consists of all other Fock states involving single- 
particle orbitals outside the mean-field space The 
inclusion of this complementary space leads to correc¬ 
tions to the mean-held energy that are of higher order in 
1 /N relative to the contribution from the primary space 
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. We shall show in the following that this fact can be 
exploited efficiently to bridge the gap between the few- 
body and thermodynamic limits. 

A. The “unit vortex” 

We begin with the relatively simple case of L/N = 1, 
where (as described above) within the mean-field approx¬ 
imation only the m = 1 orbital is occupied. As a result 
there is a single vortex located at the centre of the trap. 
The primary space thus consists of only one Fock state, 
with all atoms occupying the m = 1 orbital. The comple¬ 
mentary space is spanned by all other Fock states that 
have a non-zero occupancy of any single-particle state 
with m 1. If one works within a truncated space in¬ 
cluding only the orbitals with to = 0,1 and 2, the yrast 
state known analytically to leading order in N 

= (2) 
k V 2 

where the ket on the right denotes the Fock state with 
single-particle states in the lowest Landau level with 
TO = 0,1 and 2, and corresponding occupancies noted 
by the exponents. Returning to the separation of the 
full Hilbert space which we described above, the primary 
space consists of the single state with k = 0, with a proba¬ 
bility 1/2, while all the other states with fc ^ 0 constitute 
the complementary space, with a probability also equal 
to 1/2. Note that the amplitudes in Eq. (I2|) decrease ex¬ 
ponentially with k, i.e., with the occupancy of the states 
belonging to the complementary space. As we will see 
below, this is a more general feature that also appears 
for larger values of the angular momentum. 

We now consider a more systematic analysis of the con¬ 
vergence of the exact solution as a function of a grad¬ 
ual increase of the contribution from the complementary 
space. In order to capture fully the finite-size correc¬ 
tions, we extend the truncated space used in Eq. ([2|) by 
including the states with 0 < to < 6 which necessitates 
the use of numerical methods. We evaluate the many- 
body states for a fixed number of particles ni in the 
TO = 1 state (which is the occupation of the primary 
space) with the remaining {N — ni) particles in the com¬ 
plementary space. A Hamiltonian matrix is constructed 
for each value of (iV — ni), and the single eigenstate of 
lowest interaction energy is selected from each matrix. A 
truncated Hamiltonian is built and diagonalized in this 
new basis of lowest-energy states to obtain the approxi¬ 
mate energy spectrum and the corresponding eigen¬ 
functions I'kX) (where the index i = 0,1,2,... labels the 
excited states), here for the example oi N = 100 atoms. 

Figure 2 shows the low-energy spectrum as evaluated 
within this scheme, as a function of the highest occu¬ 
pancy of the complementary space Nc (red circles). (The 
energy is given in units of vq = uq f |</oo|^ where (/oo 
is the single-particle ground state of the two-dimensional 
harmonic oscillator). The right side of this graph also 


shows the energy spectrum evaluated within the usual 
full diagonalization of the many-body Hamiltonian (blue 
circles) with the same truncation, 0 < to < 6. 

We see that there is a rapid convergence of the ap¬ 
proximate solution to the exact result. Low-lying excited 
states are also reproduced fairly well, with larger devia¬ 
tions being apparent only in the higher-energy section of 
the spectrum. The relative error between the eigenener- 
gies as evaluated within our model and with full diago¬ 
nalization (in the same subspace) decreases exponentially 
with Nc- Remarkably, our method reproduces the yrast 
state as well as the low-lying excited states with high 
accuracy. As shown in panel (a) of Fig. 3 for the yrast 
state (“Y”) and the first non-trivial excited state “G” 
that is related to the Goldstone mode [s^, relative er¬ 
rors as small as 10“^ to 10“are obtained for values 
of Nc Ri N/2. (The first excited state, labeled “GM”, 
is a trivial center-of-mass excitation). Panel (b) of Fig. 3 
shows a logarithmic plot of the deviation from unity over¬ 
lap of the model yrast state and the exact yrast state, 
1— I ('k^ I Tex) P, as a function of Nc- This plot clearly 
shows an exponential convergence with the number of 
particles in the complementary space. Finally, Fig. 3 (c) 
shows the size of the submatrices arising in our calcula¬ 
tions. Note that all of these submatrices are dramatically 
smaller than that of the full Cl matrix, which has a di¬ 
mension of 189, 509 for the specific example considered 
here. 


B. States of discrete rotational symmetry 

While the above example of the unit vortex at L/N = 1 
is instructive, it is special in the sense that within the 
mean-field approximation there is only one state that is 
macroscopically occupied. We therefore consider, with¬ 
out loss of generality, the ratio L/N = 9/5 where within 
the mean-field approximation two vortices have nucleated 
in the cloud. 

In what follows below we show how our method is 
applied considering this case as an example. We start 
with some approximate and semi-analytic results, which 
demonstrate the use of our method. At the end of this 
section we present the full numerical results. 

In order to simplify the discussion, we first truncate 
the space to the single-particle states with 0 < to < 4. 
(It is straightforward to generalize the arguments pre¬ 
sented below to larger spaces). For large N, there is a 
macroscopic occupancy of the single-particle states with 
TO = 0, 2 and 4 in the primary space. In this limit, it is 
convenient to approximate the exact yrast state by the 
sum 

|To) = ^(-l)'=^fe|0”‘>W,2”=W^4”'‘('=)), (3) 

k 

where the occupancy of the orbitals with to = 0, 2 and 
4 are n^ik), n 2 (k) and n 4 {k), respectively. Fixing the 
occupancy of one orbital, for example, no{k) = k, the 
constraints of fixed particle number and fixed total an¬ 
gular momentum determine the occupancy of the other 
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FIG. 2: Interaction energies (in units of vo, see text) for 
L = N = 100 evaluated within the model (left, red circles) as 
a function of the highest occupancy Nc of the complementary 
space for a single-particle basis of orbitals with 0 < m < 6. 
The blue bullets to the right show the result from the “ex¬ 
act” approach in the same subspace of single-particle states 
(see text). The relative error between the model and the full 
diagonalization in the same subspace is given in panel (b) 
of Fig. 3 for the yrast state (“Y”) and the Goldstone mode 
(“G”) (see text), converging to numbers as small as 10“^ to 
IQ-io j-Qj. vaj^es of Nc « N/2. (The state labeled ”GM” is a 
center-of-mass excitation). 


two. Consequently, for the above choice of nQ^k), we have 
n 2 {k) = llA^/10 — 2k and n 4 (fc) = k — N/10. As dis¬ 
cussed in Ref. , within the states contained in Eq. Q 
the eigenvalue equation takes the form 

— + V'fc.feV'fe ~ Vk^k+ii’k+i = E-ipk- (4) 

Here 14,fe' are the matrix elements of the interaction V 
between the states |fc) and |fc') on the right of Eq. ([3|), 
and E is the interaction energy. The above equation 
is obvious, since in addition to the diagonal matrix ele¬ 
ments, the interaction, being a two-body operator, con¬ 
nects states where two atoms from the states with m = 0 
and m = 4 get transferred to the state with m = 2, and 
vice versa. 

Assuming that is a smooth and differentiable func¬ 
tion of fc, this eigenvalue equation assumes the familiar 
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FIG. 3: (a) Relative error in the interaction energy, and (b) 
the deviation from unity overlap of the model yrast state with 
the exact yrast state, 1— | (^a I ^ex) |^, on a logarithmic scale 
as a function of the highest occupancy of the complementary 
space, Nc, for the data of Fig. 2. Panel (c) shows the dimen¬ 
sion of the submatrices as a function of Nc. 


form 


- + Vk^k+i)dltlJk + Ves{k)'ipk = Eipk, (5) 

where Ves{k) = Vk,k - 14,fe-i - 14,fe-i-i- In the vicinity of 
its minimum at some k^ where Veff has the value Eq, this 
effective potential can be approximated as Veff — Eq (x 
{k — ko)^/2. Thus, 'i/'fe satisfies an eigenvalue equation of 
a linear single-particle problem in the effective potential 
Ves{k). Eor the example of L/N = 9/5, in the limit of 
very large N we find that the minimum of Ves{k) occurs 
at ko/N = uq/N Ri 0.2289, a number that is very close to 
the value of the mean-field coefficient for |coP(« 0.2307) 
given in the Table above (albeit here for 0 < m < 4). 
The energy of the yrast state within the primary space is 
found to be A R! (0.18807V^ — 0.31497V)tio plus terms of 
order unity. This expression for E must now be corrected 
at 0{N) for the effects of the complementary space. 

To gain insight into the role of the orbitals outside the 
mean-field space, we first consider contributions to the 
complementary space due to the m = 1 single-particle 
state only. The ipk in Eq. ([3]) for the yrast state is a 
Gaussian with a width of 0{'\/N). If one promotes 2ni 
particles to the single-particle state with m = 1, where 
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ni is of 0{N^), to leading order in N, the corresponding 
yrast state can be written as 

| 2 ni) oc (a|aiaoa 2 )"" I'f'o). ( 6 ) 

Here |2ni) denotes the yrast state with 2ni atoms in the 
single-particle state with m = 1 and aj„ and am are the 
usual creation and annihilation operators for a particle 
with angular momentum mTi. The Gaussian form of the 
amplitudes of the primary space components of the state 
| 2 ni) is preserved, and its center is simply shifted by a 
term of order unity. As before, the primary components 
are all of 0{N) with a width of 0{'\fN), and the oc¬ 
cupancy of the states in the complementary space is of 
0{N^). As a consequence, the energy of the state |2ni) 
is the same as that of |'I'o) to leading order in N, i.e., to 
0{N^), and it is only necessary to consider corrections 
from the complementary space which are of subleading 
order, i.e., 0{N). This implies that it is sufficient to 
approximate the full state of Eq. (|31) by the single com¬ 
ponent l^-o) = This is a very 

considerable simplification. 

Using this single component (appropriately renormal¬ 
ized to unity), we find that neither the curvature of 14 ff 
nor the “inertial parameter” (Vk^k-i + 14 ,fc-i-i) depend 
on ni to leading order. Further, we find that the diago¬ 
nal energies scale linearly with ni, while the off-diagonal 
matrix elements (2ni, 0|U|2ni -I- 2, 0) are seen to be pro¬ 
portional to N (2ni -|- l)(2ni - 1 - 2). These are the only 
nonzero off-diagonal matrix elements which come from 
the operator aioiajaj (plus its Hermitian conjugate). 
Diagonalizing the resulting matrix, we find that the yrast 
energy becomes E « (0.18807V^ — 0.3885A^)r;o and that 
the probabilities of the various states with 2 ni atoms in 
the m = 1 state decrease exponentially with ni. We 
emphasize here that the above result is generic and not 
specific for the example that we have considered (see Ap¬ 
pendix). 

Let us now turn to the case where contributions from 
both the m = 1 and to = 3 single-particle orbitals are in¬ 
cluded in the complementary space. Generalizing Eq. m, 
we see that the states including contributions from the 
complementary space can be constructed as 

| 2 ni + 1 , 2^3 + 1 ) oc l^-o), 

| 2 ni, 2 n 3 ) oc afl^o). (7) 

We emphasize that the counting of states in this equa¬ 
tion is correct: Naively, one might expect that two dis¬ 
tinct states, e.g., with ni = = 1 could be constructed 

from |d>o) through the application of either 0 ^ 030202 , or 
a|a 3 aoa 4 . Given the Gaussian nature of tjjk, however, 
these states are not orthogonal, and only one of them 
should be included. 

There are four classes of nonzero off-diagonal matrix 
elements. The first two classes include 

(ni -I- 2,7131 U ^^ 3 ) = ?^2Uii,2o y /(ni -|- l)(ni -|- 2), 


and 

(rii, 773 -I- 2|U|ni, 773 ) = y/n^rilV'iz,24: yj (773 -I- 1)(773 -I- 2), 

where Vij^ki = {(t>Qi,4>oj\V\(t>ok,(l>m)- The remaining two 
classes reflect the staggering of the ground state wave 
function (i.e., the alternating sign seen explicitly in Eq.[2] 
above). They have the form 

(77i, 773|U|77i - 1,773 -f 1) = (“ 1)”^ \/77l (773 -|- 1) 

X (2y^77o 772Ui 2,30 — 2y772T74Ui4_3o) (8) 

and 

(771 -h 1,773 -h ^V\ni,ny) = (-l)"^v^(77i -h 1)(773 -h 1) 

X (772 V 22 .13 — 2-^770 7741/04, 13 )- (9) 

For Nc particles in the complementary space, the ma¬ 
trix has dimension Nc/2 + \. Gonsidering as an example 
Nc = 18, the yrast energy E Ri (0.1880iV^ — O.SldQTvjuo 
obtained above for the primary space is shifted by the 
amount AE = —0.1893Auo- As previously discussed for 
the simpler case with only the to = 1 orbital from the 
complementary space, the amplitudes show an exponen¬ 
tial decrease with 774 and 773 . 

In the preceding paragraphs we approximated the 
yrast state with the single component [To) = 
|Qno(feo)^ 4”'*^^“)) which allowed us to give a de¬ 

tailed description of our method, providing us with some 
semi-analytic results. In what follows below we apply 
the same method fully numerically, without making the 
approximation for the yrast state consisting of a single 
component and compare it with the result of the full diag- 
onalization (in the spirit of the case of L/N = 1 discussed 
above). Here, however, we consider only the yrast state, 
since the analysis of the excitation spectra beyond the 
value of L/N = 1 is more complicated due to the fact 
that, even for the low-lying states, competing solutions 
of different symmetry are found, see Ref. (43 |. 

We follow the same procedure as above, consider¬ 
ing two different truncations of the Hilbert space with 
0 < TO < 4 and 0 < to < 6. First, we evaluate the many- 
body states of Eq. 0 for each configuration of particles 
in the complementary space. Then, the lowest-energy 
state for each such configuration is chosen, and the trun¬ 
cated Hamiltonian is diagonalized in this new basis of 
lowest-energy states to obtain the approximate ground 
state energy E\ and corresponding wavefunction 
Figure 4 shows the quantity 1 — [Tex)!^ versus N in 

the space 0 < to < 4 (squares) and 0 < to < 6 (circles), 
on a double-logarithmic scale. Intriguingly, for large 
N, this quantity shows a simple power-law behaviour, 
l-|(^°4|Tex)PcX A-U 

The fact that the overlap becomes worse for the more 
extended space, 0 < to < 6 , as compared to the more 
restricted one, 0 < to < 4, is not a surprise, since the rel¬ 
ative dimensionality of the full Hilbert space compared 
with the dimensionality of the space of our approach in¬ 
creases dramatically in the more extended space. The 
crucial result of this analysis is that the Fock state am¬ 
plitudes for a given configuration in the complementary 
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FIG. 4: The quantity 1 — as a function of N for 

a fixed angular momentum per particle, L/N = 9/5. Here, at 
each point the number Nc of particles in the complementary 
space was increased until numerical convergence of the result 
was obtained. The red squares refer to the space with 0 < 
m < 4 (open squares, K = 1, and solid squares, JsT = 4, 
where K is the number of lowest-energy states included in 
the reduced matrix), and the blue circles to the space with 
0 < m < 6 (open circles, K = 1, and solid circles, isT = 4). 


space decrease exponentially with the number of particles 
in that space in both cases (we also emphasize that the 
number of the single-particle orbitals which are occupied 
by a macroscopic number of atoms saturates quickly for 
larger values of m, as it is clearly seen from the logarith¬ 
mic plot of the occupancies in Fig. 5). The exponent is 
independent of N in the limit N ^ oo, indicating that 
the number of particles in the complementary space has 
a limit of 0{N^). This behavior is in fact strongly sup¬ 
ported by the occupancies obtained through direct diag- 
onalization, as seen also in Fig. 5. The occupancy of the 
states which construct the primary space in compliance 
with the mean-field approximation (red lines) increases 
(except for the orbital with m = 6 that saturates), while 
the orbitals of the complementary states show a power- 
law decay with the number of particles. 

More generally, one can retain the K eigenstates of 
lowest energy for each configuration of particles in the 
complementary space. (This corresponds to retaining ex¬ 
cited states of the effective harmonic-oscillator problem 
in the primary space). For a given basis set, an increase 



Particle Number N 

FIG. 5: Occupancies of the exact yrast states at L/N = 9/5 as 
a function of the particle number, N, on a double-logarithmic 
scale here for a basis with 0 < m < 7. The saturation of the 
single-particle basis is reflected by the significant reduction 
in occupancy for higher values of m. The red circles indicate 
single-particle states that belong to the primary space; the 
green squares indicate those that belong to the complemen¬ 
tary space. Clearly, with increasing N there is a convergence 
towards the occupancies obtained from mean-field (see the 
Table above). 

of K accelerates the convergence towards the full solution 
(see Fig. a. For fixed N, one finds exponential conver¬ 
gence in K. In this case, however, the exponent depends 
on N such that convergence is more rapid for larger par¬ 
ticle numbers, and the generalized oscillator ground state 
alone contributes to the yrast state as iV —>■ oo. We stress 
that for Nc = N and retaining all the possible K states, 
the approach is simply a passive unitary transformation 
of the basis and the results are necessarily identical to 
the full exact solution. 


IV. CONCLUSIONS 

In short, this paper suggests a significantly simpli¬ 
fied understanding of the properties of a rotating Bose- 
Einstein condensate of trapped atoms. The direct nu¬ 
merical strategy for this problem would be to include a 
certain set of single-particle states and to diagonalize the 
resulting many-body Hamiltonian matrix. The difficulty 
is that the dimension of this matrix grows prohibitively 
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as the number of particles or the angular momentum 
(and thus, consequently, the number of necessary single¬ 
particle basis states) increases. 

The method presented in this study makes use of the 
fact that only certain single-particle states are macro- 
scopically occupied, while all other states have an occu¬ 
pancy of order unity. This introduces a natural sepa¬ 
ration of the Hilbert space into a primary and a com¬ 
plementary part. The first, containing the macroscopi- 
cally occupied single-particle states, can be regarded as 
a generalised harmonic oscillator problem that can with¬ 
stand major truncation when N is large. The resulting 
simplification is significant: The size of the Hamiltonian 
matrix can be reduced safely by a factor of order 
where k is the number of single-particle states included in 
the primary space. The contribution of the complemen¬ 
tary subspace to the many-body states falls exponentially 
with the number of particles in it. Therefore, the present 
approach shows clearly that the vast majority of these 
states do not make a significant contribution to the yrast 
states, providing a simple understanding of how scatter¬ 
ing processes between the primary and complementary 
spaces govern the transition from finite-sizes to the ther¬ 
modynamic limit. 

At mean-field level there are discontinuous phase tran¬ 
sitions between states of different rotational symmetry 
corresponding to level crossings . The states involved 
in such crossings can be constructed using the meth¬ 
ods described here for distinct choices of the primary 
space. In the immediate vicinity of the crossing point, 
these states will be nearly degenerate and can in prin¬ 
ciple mix. Fortunately, however, the fact that they are 
based on different primary spaces ensures that matrix el¬ 
ements of the interaction between these states will vanish 
exponentially with the number of particles. Thus, such 
mixing becomes increasingly unimportant as the number 
of particles grows unless one is precisely at the crossing 
point. 

The analysis presented here has been restricted to the 
limit of weak interactions, where one may neglect the 
single-particle eigenstates of the harmonic-oscillator po¬ 
tential with radial nodes. We stress, however, that our 
results are quite general and are not specific to this per¬ 
turbative regime. Rather, the illustration of our method 
for the case of weak interactions represents a “proof of 
principle” and provides a representative example of our 
approach. Even in the regime of stronger interactions, 
which is actually of greater experimental relevance, one 
can solve the mean-field Gross-Pitaevskii equation to de¬ 
termine which states should be contained in our primary 
subspace. This subspace could, in principle, contain any 
or even all Landau levels with an angular momentum 
consistent with the discrete rotational symmetry of the 
ground state. In practice, mean-field calculations for 
stronger interactions show that the probability of find¬ 
ing states with Ur radial nodes in the mean-field wave 
function decreases exponentially with increasing Ur- The 
inclusion of additional Landau levels will not lead to any 
material complication in the construction or the solution 
of the generalized harmonic-oscillator problem described 


in Sec. HI A. Excited Landau levels for other angular 
momenta will contribute to the complementary space. 
Since the exponential convergence found in the present 
manuscript is dictated by angular momentum considera¬ 
tions and not by the radial structure of the single-particle 
wave functions, this feature will be unaltered by the in¬ 
clusion of higher Landau levels. In short, the scheme 
introduced here will remain valid and useful even if addi¬ 
tional Landau levels are included. However, as mentioned 
above, internal convergence criteria must be adopted for 
assessing the accuracy of such calculations since full nu¬ 
merical diagonalizations will not be practical. 

While the truncations adopted here appear to be par¬ 
ticularly promising when the number of particles is large, 
we have not proven that the approach is a viable quan¬ 
titative alternative to the exact diagonalization for large 
systems. Such a claim would require extensive bench¬ 
marking that is beyond the scope of this work and re¬ 
mains a matter for further investigation. Rather, the 
present study offers new and very explicit insight into the 
structure of the many-body wavefunction and its relation 
to the mean-field approximation for a rotating atomic 
superfluid. Although we have focused primarily on the 
yrast line, the procedure adopted here should also be 
suitable for investigating the richness of the excitation 
spectrum. We stress that the method is physically well- 
motivated and provides a well-defined transformation of 
the basis of many-body states that is completely passive 
in the sense that it suggests the order but not the degree 
of truncations of the basis. 
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APPENDIX: Toy Model 

Consider, e.g., a real symmetric matrix which is zero 
except for the diagonal matrix elements = n — 1 and 
the off-diagonal matrix elements with An^n+i = 

nf. The eigenvalue equation 

An,n—l(-n—l “f ^n,nCn — ECji 

has the form 

(n - l)/c„_i + {n- l)c„ -I- nf Cn+i = Ec„. 

The lowest eigenvalue of this matrix has a finite value in 
the limit that its dimension approaches infinity, provided 
that 0 < I/I < 1/2. Subject to this restriction, the lowest 
eigenvalue is E = — |/|a: and the corresponding solution 
is Cn = x^\/l — with X = —(1 — \/l — Ap)/{2f). For 
large values of n, the above eigenvalue equation has the 
simpler form fcn-i + Cn + fcn+i = 0, which has the 




solution c„ oc cc", independent of i?. In other words, c„ exponential convergence, 
decays exponentially at precisely the same rate for all 
eigenvectors of finite energy. Although elementary, this 
toy model illustrates the present mechanism leading to 
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